{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Global Uncertainty Analysis: Polynomial Chaos Expansion (PCE) for Chemical Reaction Systems\n",
    "\n",
    "\n",
    "This ipython notebook uses MUQ as a basis for adaptive Polynomial Chaos Expansions to perform global uncertainty analysis for chemical reaction systems.  This ipython notebook details a workflow using RMG, Cantera, and MUQ codes.\n",
    "\n",
    "Muq binary only works on linux systems, please also add the ~/anaconda/envs/your_env/lib folder to your $PYTHONPATH to import muq smoothly."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "from rmgpy.tools.canteraModel import Cantera, getRMGSpeciesFromUserSpecies\n",
    "from rmgpy.species import Species\n",
    "from rmgpy.chemkin import loadChemkinFile\n",
    "from rmgpy.tools.muq import ReactorPCEFactory\n",
    "from rmgpy.tools.uncertainty import Uncertainty"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Load chemical kinetic mechanism from RMG chemkin file and dictionary."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# The paths for the chemkin file and its species dictionary\n",
    "chemkinFile = 'uncertainty/chem_annotated.inp'\n",
    "dictFile = 'uncertainty/species_dictionary.txt'"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For uncorrelated uncertainty studies, we can simply create a cantera job for the model using species and reactions lists from the `loadChemkinFile` function.  Alternatively the `Uncertainty` class's `loadModel` function can be used."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "outputDir = 'uncertaintyUncorrelated'\n",
    "speciesList, reactionList = loadChemkinFile(chemkinFile, dictFile)\n",
    "\n",
    "# Declare some species that we want to use as initial conditions or in our uncertainty analysis\n",
    "PDD = Species().fromSMILES(\"CCCCCCCCCCCCc1ccccc1\")\n",
    "C11ene=Species().fromSMILES(\"CCCCCCCCCC=C\")\n",
    "ETHBENZ=Species().fromSMILES(\"CCc1ccccc1\")\n",
    "\n",
    "# Map the species to their respective objects in the speciesList and reactionList\n",
    "mapping = getRMGSpeciesFromUserSpecies([PDD,C11ene,ETHBENZ], speciesList)\n",
    "\n",
    "reactorTypeList = ['IdealGasConstPressureTemperatureReactor']\n",
    "molFracList = [{mapping[PDD]: 1.0}]\n",
    "Tlist = ([623],'K')\n",
    "Plist = ([350],'bar')\n",
    "reactionTimeList = ([72], 'h')\n",
    "\n",
    "# Create the cantera model\n",
    "job = Cantera(speciesList=speciesList, reactionList=reactionList, outputDirectory=outputDir)\n",
    "# Load the cantera model based on the RMG reactions and species\n",
    "job.loadModel()\n",
    "# Generate the conditions based on the settings we declared earlier\n",
    "job.generateConditions(reactorTypeList, reactionTimeList, molFracList, Tlist, Plist)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's create a second cantera model and output folder to use for the correlated uncertainty analysis.  Be careful to use NEW species and reaction objects, since running the two cantera models through the global uncertainty analysis module at the same time inside a single python script may cause collision issues (the ReactorPCEFactory class actively manipulates species and reaction object data).  For an uncorrelated analysis, we must use the `Uncertainty` class, because we need to pass in the partial input uncertainty dictionaries: `Uncertainty.kineticInputUncertainties` and `Uncertainty.thermoInputUncertinaties` into the global analysis uncertainty class."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "outputDir = 'uncertaintyCorrelated'\n",
    "\n",
    "uncertainty = Uncertainty(outputDirectory='testUncertainty')\n",
    "uncertainty.loadModel(chemkinFile, dictFile)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Use the Uncertainty class to assign correlated uncertainties, which we will propagate.  This requires an additional step of loading the RMG database and extracting the original parameter sources (i.e. rate rules and thermo)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "uncertainty.loadDatabase()\n",
    "uncertainty.extractSourcesFromModel()\n",
    "# Assign correlated parameter uncertainties \n",
    "uncertainty.assignParameterUncertainties(correlated=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Again create a Cantera model object that stores the reaction conditions."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "mappingCorrelated = getRMGSpeciesFromUserSpecies([PDD,C11ene,ETHBENZ], uncertainty.speciesList)\n",
    "\n",
    "reactorTypeList = ['IdealGasConstPressureTemperatureReactor']\n",
    "molFracList = [{mappingCorrelated[PDD]: 1.0}]\n",
    "Tlist = ([623],'K')\n",
    "Plist = ([350],'bar')\n",
    "reactionTimeList = ([72], 'h')\n",
    "\n",
    "jobCorrelated = Cantera(speciesList=uncertainty.speciesList, reactionList=uncertainty.reactionList, outputDirectory=outputDir)\n",
    "# Load the cantera model based on the RMG reactions and species\n",
    "jobCorrelated.loadModel()\n",
    "# Generate the conditions based on the settings we declared earlier\n",
    "jobCorrelated.generateConditions(reactorTypeList, reactionTimeList, molFracList, Tlist, Plist)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Input a set of kinetic $(k)$ and thermo $(G)$ parameters to be propagated and their uncertainties $(d\\ln(k), dG)$ into the `ReactorPCEFactory` class.  These kinetic and thermo parameters should typically be pre-screened from local uncertainty analysis to narrow down to the most influential parameters.  \n",
    "Each parameter's uncertainty is considered to be a uniform uncertainty interval where unit random variables $\\ln(k)_{rv}$ and $G_{rv}$ are mapped to the user-assigned parameter uncertainties.\n",
    "\n",
    "$\\ln(k)_{rv} \\sim U(-1, 1) \\rightarrow ln(k) \\sim U(-d\\ln(k), d\\ln(k))$\n",
    "\n",
    "$G_{rv} \\sim U(-1, 1) \\rightarrow G \\sim U(-dG, dG)$\n",
    "\n",
    "\n",
    "Polynomial chaos expansions (PCE) are contructed for the desired outputs of interest (species mole fractions)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In the case of correlated parameters, we will need to indicate them in the `ReactorPCEFactory` object, and as well as provide the dictionaries of assigned partial uncertainties from the `Uncertainty` class."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# Create ReactorPCEFactory global uncertainty analysis object for the uncorrelated case\n",
    "\n",
    "reactorPCEFactory = ReactorPCEFactory(cantera=job,\n",
    "                            outputSpeciesList=[mapping[PDD], mapping[C11ene]],\n",
    "                            kParams=[796, 801], # [Styrene+Decyl=Rad1, c10ene + ebzyl = rad4]\n",
    "                                    # A list of indices corresponding to the uncertain reactions\n",
    "                            kUncertainty = [1.414, 1.414],   \n",
    "                                      # a list of dlnk corresponding to uncertainties of kParams\n",
    "                            gParams = [4,5],  # A list of indices corresponding to the uncertain thermo   [PDD, Toluene]\n",
    "                            gUncertainty = [2.0, 1.0],\n",
    "                                      # a list of values corresponding to dG\n",
    "                            correlated = False \n",
    "                            )\n",
    "\n",
    "reactorPCEFactoryCorrelated = ReactorPCEFactory(cantera=jobCorrelated,\n",
    "                            outputSpeciesList=[mappingCorrelated[PDD], mappingCorrelated[C11ene]],\n",
    "                            kParams=['R_Addition_MultipleBond Cds-HH_Cds-CsH;CsJ-CsHH', \n",
    "                                    'Estimation STYRENE(3)+DECYL(56)=RAD1(14)'],   \n",
    "                                      # labels for the correlated kinetics parameters\n",
    "                            kUncertainty = uncertainty.kineticInputUncertainties,   \n",
    "                                      # a list of dictionaries that gives the reaction's partial uncertainties\n",
    "                                      # with respect to the string labels of the kinetic correlated parameters, i.e. 'H_Abstraction CHO/Oa'\n",
    "                            gParams = ['Group(group) Cs-CsCsHH', 'Library TOLUENE(2)'],  # labels for the correlated thermo parameters\n",
    "                            gUncertainty = uncertainty.thermoInputUncertainties,\n",
    "                                      # a list of dictionaries that gives the species partial uncertainties\n",
    "                                      # with respect to the string labels of the correlated thermo parameters, i.e. 'Group(ring) cyclohexane'\n",
    "                            correlated = True   \n",
    "                            )\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Begin generating the PCEs adaptively based a runtime.\n",
    "\n",
    "There are actually three methods for generating PCEs. See the `ReactorPCEFactory.generatePCE` function for more details.\n",
    "\n",
    "- Option 1: Adaptive for a pre-specified amount of time\n",
    "- Option 2: Adaptively construct PCE to error tolerance\n",
    "- Option 3: Used a fixed order, and (optionally) adapt later.  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "reactorPCEFactory.generatePCE(runTime=60)  # runtime of 60 seconds."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's compare the outputs for a test point using the real model versus using the PCE approximation.\n",
    "Evaluate the desired output mole fractions based on a set of inputs `ins = [[` $\\ln(k)_{rv}$ `], [` $G_{rv}$ `]]` which contains the \n",
    "random unit uniform variables attributed to the uncertain kinetics and free energy parameters, respectively."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "trueTestPointOutput, pceTestPointOutput = reactorPCEFactory.compareOutput([0.5,0.2,0.1,-.7])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Obtain the results: the species mole fraction mean and variance computed from the PCE, as well as the global sensitivity indices."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "mean, variance, covariance, mainSens, totalSens = reactorPCEFactory.analyzeResults()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Do the same analysis for the correlated `reactorPCEFactory`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "reactorPCEFactoryCorrelated.generatePCE(runTime=60)  # runtime of 60 seconds."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "trueTestPointOutput, pceTestPointOutput = reactorPCEFactoryCorrelated.compareOutput([0.5,0.2,0.1,-.7])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "mean, variance, covariance, mainSens, totalSens = reactorPCEFactoryCorrelated.analyzeResults()"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 2",
   "language": "python",
   "name": "python2"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 2
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython2",
   "version": "2.7.11"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}
